

####
# Replication code for "Mobilizing from Below: The Home State Effect on Climate Coalitions"
# Authors: Jonas Meckling; Samuel Trachtman
# Governance
####

rm(list=ls())
setwd("C:/Users/sam.trachtman/Dropbox/Climate federalism/replication")

#libraries
library(lmtest)
library(tidyverse)
library(ggplot2)
library(psych)
library(rio)
library(arm)
library(dplyr)
library(sandwich)

##### GENERATING POLICY MEASURES #####

#load policy data 
load("policydat.Rdata")

#generate mean policy score over 10 year window 
fa_mean <- three_measures %>%
  group_by(abb) %>%
  dplyr::summarize(
    rps_score = mean(rps_over_tot,na.rm=T),
    pricing_score = mean(capped,na.rm=T),
    ee_score = mean(electric_spend,na.rm=T)) #%>%

#estimate PCA
fa_mean_analysis <- fa_mean %>% dplyr::select(-abb)
pca <- prcomp(fa_mean_analysis, scale = T, center = T)
pca$rotation

#pull pca scores 
fa_mean$scores <- pca$x[,1]
hist(fa_mean$scores)

#normalize from 0 to 1 
fa_mean$scores_mean <- (fa_mean$scores - min(fa_mean$scores))/(max(fa_mean$scores) - min(fa_mean$scores))
names(fa_mean) <- c("abb", "rps_mean", "pricing_mean", "ee_mean", "raw_score_mean", "scores_mean")

#plotting Figure A1 
policy.plot <- ggplot(data = fa_mean,
                      mapping = aes(x = scores_mean, y = reorder(abb, scores_mean)))+
  geom_vline(xintercept = 0, color = "gray30") +
  geom_point(size = 1.5)+
  theme(axis.title = element_blank())+
  xlab("")+
  ylab("")+
  theme_bw()
policy.plot

##### PREPPING DATA FOR ANALYSIS #####

#load utility data 
load("utility_state_level.Rdata")

#grouping by utility dataset 1 for electricity distribution weights 
agg.df.ut = df %>%
  group_by(utilityid) %>%
  summarise_at(vars(tot_control_num, state_position, pres_dem_2016, 
                    scores_2015, rps_2015, pricing_2015, ee_2015, 
                    scores_mean, rps_mean, pricing_mean, ee_mean, 
                    happening, gdp_per, per_capita_tax, wind_potential_adj, solar_potential_adj, medianaqi, unemp, electricity_price), funs(weighted.mean(., sales_mwh)))
names(agg.df.ut) <- c("utilityid", "control.ut", "st.position.ut", "dem.ut", 
                      "scores_2015.ut", "rps_2015.ut", "pricing_2015.ut", "ee_2015.ut", 
                      "scores_mean.ut", "rps_mean.ut", "pricing_mean.ut", "ee_mean.ut", 
                      "happening.ut", "gdp_per.ut", "per_capita_tax.ut", "wind_potential_adj.ut", "solar_potential_adj.ut", "medianaqi.ut", "unemp.ut", "electricity_price.ut")

#grouping by utility dataset 2 for electricity generation weights 
agg.df.gen = df %>%
  group_by(utilityid) %>%
  summarise_at(vars(tot_control_num, state_position, pres_dem_2016, 
                    scores_2015, rps_2015, pricing_2015, ee_2015, 
                    scores_mean, rps_mean, pricing_mean, ee_mean, 
                    happening, gdp_per, per_capita_tax, wind_potential_adj, solar_potential_adj, medianaqi, unemp, electricity_price), funs(weighted.mean(., capacityall)))
names(agg.df.gen) <- c("utilityid", "control.gen", "st.position.gen", "dem.gen", 
                       "scores_2015.gen", "rps_2015.gen", "pricing_2015.gen", "ee_2015.gen", 
                       "scores_mean.gen", "rps_mean.gen", "pricing_mean.gen", "ee_mean.gen", 
                       "happening.gen", "gdp_per.gen", "per_capita_tax.gen", "wind_potential_adj.gen", "solar_potential_adj.gen", "medianaqi.gen", "unemp.gen", "electricity_price.gen")

#collapse utility-level variables to parent level 
agg.df <- df %>%
  dplyr::group_by(utilityid, utility, parent) %>%
  dplyr::summarize(perc_coal = sum(capacitycoal)/sum(capacityall),
                   perc_renew = sum(capacityrenewable)/sum(capacityall),
                   capacitycoal = .001*sum(capacitycoal),
                   capacityrenewable = .001*sum(capacityrenewable),
                   capacitytot = sum(capacityall),
                   sales_mwh = sum(sales_mwh),
                   cpp_support = mean(direct_cpp_support),
                   cpp_oppose = mean(direct_cpp_oppose),
                   uarg = mean(uarg2017),
                   oppose_coalition = mean(oppose_coalition),
                   ownership = max(ownership2),
                   num_oppose = mean(num_oppose),
                   num_support = mean(num_support),
                   acore = mean(acore)
  )
agg.df$perc_coal[is.nan(agg.df$perc_coal)] <- 0 #NA to 0
agg.df$perc_renew[is.nan(agg.df$perc_renew)] <- 0 #NA to 0 

#merge for the utility-level and state-level independent variables 
agg.df <- left_join(agg.df, agg.df.gen, by = "utilityid")
agg.df <- left_join(agg.df, agg.df.ut, by = c("utilityid"))

#coding CPP outcome as factor 
agg.df$cpp <- factor(ifelse(agg.df$cpp_support==1, "support",
                            ifelse(agg.df$cpp_oppose==1, "anti", "neutral")))

#cpp alternative CPP outcome as factor 
agg.df$cppc <- factor(ifelse(agg.df$cpp_support==1, "support",
                             ifelse(agg.df$oppose_coalition==1, "anti",
                                    ifelse(agg.df$cpp_oppose==1, "anti", "neutral"))))

#cpp numeric coding 
agg.df$cppc.num <- ifelse(agg.df$cpp_support==1, 1,
                          ifelse(agg.df$oppose_coalition==1, -1,
                                 ifelse(agg.df$cpp_oppose==1, -1, 0)))

#coding ownership and type of utility, types: 0 = IPP, 1 = IOU, 2 = Muni, 3 = Coop
agg.df$type.c <- with(agg.df,
                      ifelse(ownership==0, "IPP",
                             ifelse(ownership==1, "IOU",
                                    ifelse(ownership==2, "MUNI", "COOP"))))
agg.df$coop = agg.df$type.c == "COOP"
agg.df$muni = agg.df$type.c == "MUNI"

#indicator variables for electric service / electricity generation / both 
agg.df$just_sales <- agg.df$capacitytot==0
agg.df$just_gen <- agg.df$sales_mwh==0

#variable for proportion of activity that is generation 
agg.df$prop_gen <- agg.df$capacitytot*8760/(agg.df$capacitytot*8760 + agg.df$sales_mwh)

#coding treatment and covariates based on generation / distribution weights
agg.df$scores_2015 <- with(agg.df,
                           ifelse(just_gen==1, scores_2015.gen,
                                  ifelse(just_sales==1, scores_2015.ut, 
                                         prop_gen*scores_2015.gen + (1-prop_gen)*scores_2015.ut)))

agg.df$scores_mean <- with(agg.df,
                           ifelse(just_gen==1, scores_mean.gen,
                                  ifelse(just_sales==1, scores_mean.ut, 
                                         prop_gen*scores_mean.gen + (1-prop_gen)*scores_mean.ut)))

agg.df$dem <- with(agg.df,
                   ifelse(just_gen==1, dem.gen,
                          ifelse(just_sales==1, dem.ut, 
                                 prop_gen*dem.gen + (1-prop_gen)*dem.ut)))
agg.df$control <- with(agg.df,
                       ifelse(just_gen==1, control.gen,
                              ifelse(just_sales==1, control.ut, 
                                     prop_gen*control.gen + (1-prop_gen)*control.ut)))

agg.df$st.position <- with(agg.df,
                           ifelse(just_gen==1, st.position.gen,
                                  ifelse(just_sales==1, st.position.ut, 
                                         prop_gen*st.position.gen + (1-prop_gen)*st.position.ut)))

agg.df$rps_2015 <- with(agg.df,
                        ifelse(just_gen==1, rps_2015.gen,
                               ifelse(just_sales==1, rps_2015.ut, 
                                      prop_gen*rps_2015.gen + (1-prop_gen)*rps_2015.ut)))

agg.df$pricing_2015 <- with(agg.df,
                            ifelse(just_gen==1, pricing_2015.gen,
                                   ifelse(just_sales==1, pricing_2015.ut, 
                                          prop_gen*pricing_2015.gen + (1-prop_gen)*pricing_2015.ut)))

agg.df$ee_2015 <- with(agg.df,
                       ifelse(just_gen==1, ee_2015.gen,
                              ifelse(just_sales==1, ee_2015.ut, 
                                     prop_gen*ee_2015.gen + (1-prop_gen)*ee_2015.ut)))

agg.df$rps_mean <- with(agg.df,
                        ifelse(just_gen==1, rps_mean.gen,
                               ifelse(just_sales==1, rps_mean.ut, 
                                      prop_gen*rps_mean.gen + (1-prop_gen)*rps_mean.ut)))

agg.df$pricing_mean <- with(agg.df,
                            ifelse(just_gen==1, pricing_mean.gen,
                                   ifelse(just_sales==1, pricing_mean.ut, 
                                          prop_gen*pricing_mean.gen + (1-prop_gen)*pricing_mean.ut)))

agg.df$ee_mean <- with(agg.df,
                       ifelse(just_gen==1, ee_mean.gen,
                              ifelse(just_sales==1, ee_mean.ut, 
                                     prop_gen*ee_mean.gen + (1-prop_gen)*ee_mean.ut)))

agg.df$happening <- with(agg.df,
                         ifelse(just_gen==1, happening.gen,
                                ifelse(just_sales==1, happening.ut, 
                                       prop_gen*happening.gen + (1-prop_gen)*happening.ut)))

agg.df$gdp_per <- with(agg.df,
                       ifelse(just_gen==1, gdp_per.gen,
                              ifelse(just_sales==1, gdp_per.ut, 
                                     prop_gen*gdp_per.gen + (1-prop_gen)*gdp_per.ut)))

agg.df$wind_potential_adj <- with(agg.df,
                                  ifelse(just_gen==1, wind_potential_adj.gen,
                                         ifelse(just_sales==1, wind_potential_adj.ut, 
                                                prop_gen*wind_potential_adj.gen + (1-prop_gen)*wind_potential_adj.ut)))

agg.df$solar_potential_adj <- with(agg.df,
                                   ifelse(just_gen==1, solar_potential_adj.gen,
                                          ifelse(just_sales==1, solar_potential_adj.ut, 
                                                 prop_gen*solar_potential_adj.gen + (1-prop_gen)*solar_potential_adj.ut)))

agg.df$medianaqi <- with(agg.df,
                         ifelse(just_gen==1, medianaqi.gen,
                                ifelse(just_sales==1, medianaqi.ut, 
                                       prop_gen*medianaqi.gen + (1-prop_gen)*medianaqi.ut)))

agg.df$unemp <- with(agg.df,
                     ifelse(just_gen==1, unemp.gen,
                            ifelse(just_sales==1, unemp.ut, 
                                   prop_gen*unemp.gen + (1-prop_gen)*unemp.ut)))

agg.df$electricity_price <- with(agg.df,
                                 ifelse(just_gen==1, electricity_price.gen,
                                        ifelse(just_sales==1, electricity_price.ut, 
                                               prop_gen*electricity_price.gen + (1-prop_gen)*electricity_price.ut)))

agg.df$per_capita_tax <- with(agg.df,
                              ifelse(just_gen==1, per_capita_tax.gen,
                                     ifelse(just_sales==1, per_capita_tax.ut, 
                                            prop_gen*per_capita_tax.gen + (1-prop_gen)*per_capita_tax.ut)))


#summary statistics 
agg.df$log_cap <- log(agg.df$capacitytot+1)
agg.df$log_sales <- log(agg.df$sales_mwh+1)
agg.df$gdp_per2 <- agg.df$gdp_per/1000
agg.df$per_capita_tax2 <- agg.df$per_capita_tax/1000

sum.stats <- data.frame(subset(agg.df, select = c(scores_mean , dem , st.position, control, perc_coal , perc_renew , log_cap , log_sales , muni, 
                                                  gdp_per2, per_capita_tax2, wind_potential_adj, medianaqi, unemp, electricity_price)))

##### REGRESSION ANALYSES FOR CPP ##### 

#descriptive relationship 
agg.df$cpp_oppose2 <- agg.df$cpp_oppose==1 | agg.df$oppose_coalition==1
mean(agg.df$cpp_oppose[agg.df$scores_mean < median(agg.df$scores_mean,na.rm=T)],na.rm=T)
mean(agg.df$cpp_support[agg.df$scores_mean < median(agg.df$scores_mean,na.rm=T)],na.rm=T)
mean(agg.df$cpp_oppose[agg.df$scores_mean > median(agg.df$scores_mean,na.rm=T)],na.rm=T)
mean(agg.df$cpp_support[agg.df$scores_mean > median(agg.df$scores_mean,na.rm=T)],na.rm=T)

#Table 1 regressions 
ord1 <- polr(cpp ~ scores_mean,
             data = agg.df)
summary(ord1)
coeftest(ord1, vcov=vcovCL(ord1, factor(agg.df$parent)))
exp(coef(ord1))

ord2 <- polr(cpp ~ scores_mean + dem + st.position  ,
             data = agg.df)
summary(ord2)
exp(coef(ord2))
coeftest(ord2, vcov=vcovCL(ord2, factor(agg.df$parent)))

agg.df$scores2 <- agg.df$scores_mean*10

ord3 <- polr(cpp ~ scores_mean +  dem + st.position + perc_coal + perc_renew + log_cap + log_sales + muni,
             data = agg.df)
summary(ord3)
coeftest(ord3, vcov=vcovCL(ord3, factor(agg.df$parent)))
exp(coef(ord3))

ord4 <- polr(cpp ~ scores_mean + dem + st.position + perc_coal + perc_renew + log_cap + log_sales + muni  
             + gdp_per2 + per_capita_tax2 + wind_potential_adj + medianaqi + unemp + electricity_price, #excluded solar potential due to missingness. 
             data = agg.df)
summary(ord4)
nobs(ord4)
coeftest(ord4, vcov=vcovCL(ord4, factor(agg.df$parent) ))



##### REGRESSION ANALYSIS FOR AD HOC COALITIONS #####

#recoding outcome
agg.df$support <- agg.df$num_support > 0 & agg.df$num_oppose <= 0
agg.df$neutral <- agg.df$num_support > 0 & agg.df$num_oppose > 0
agg.df$oppose2 <- agg.df$num_oppose > 0
agg.df$coalition_cat <- with(agg.df,
                             ifelse(num_support > 0 & num_oppose == 0, 1,
                                    ifelse(num_support==0 & num_oppose > 0, -1, 0)))
agg.df$coalition_fac <- as.factor(agg.df$coalition_cat)

#descriptive relationship 
agg.df$cpp_oppose2 <- agg.df$cpp_oppose==1 | agg.df$oppose_coalition==1
mean(agg.df$cpp_oppose2[agg.df$scores_mean < median(agg.df$scores_mean,na.rm=T)],na.rm=T)
mean(agg.df$cpp_support[agg.df$scores_mean < median(agg.df$scores_mean,na.rm=T)],na.rm=T)
mean(agg.df$cpp_oppose2[agg.df$scores_mean > median(agg.df$scores_mean,na.rm=T)],na.rm=T)
mean(agg.df$cpp_support[agg.df$scores_mean > median(agg.df$scores_mean,na.rm=T)],na.rm=T)

#table 2 regressions 
agg.df$dem2 <- agg.df$dem*10 #recode for interpretability 

reg.1 <- polr(coalition_fac ~ scores_mean,
              data = agg.df)
summary(reg.1)
coeftest(reg.1, vcov=vcovCL(reg.1, factor(agg.df$parent)))


reg.2 <- polr(coalition_fac ~ scores_mean + dem + control,
              data = agg.df)
coeftest(reg.2, vcov=vcovCL(reg.2, factor(agg.df$parent)))

reg.3 <- polr(coalition_fac ~ scores_mean + dem + control + 
                perc_coal + perc_renew + log_cap + log_sales + muni,
              data = agg.df)
coeftest(reg.3, vcov=vcovCL(reg.3, factor(agg.df$parent)))

reg.4 <- polr(coalition_fac ~ scores_mean + dem + control + perc_coal + perc_renew + log_cap + log_sales + muni 
              + gdp_per2 + per_capita_tax2 + wind_potential_adj + medianaqi + unemp + electricity_price,
              data = agg.df)
summary(reg.4)


#table A5 regressions
agg.df$member_both <- agg.df$num_support > 0 & agg.df$num_oppose > 0

reg.1 <- polr(coalition_fac ~ scores_mean,
              data = subset(agg.df,member_both==0))
summary(reg.1)

reg.2 <- polr(coalition_fac ~ scores_mean + dem + control,
              data = subset(agg.df,member_both==0))
summary(reg.2)

reg.3 <- polr(coalition_fac ~ scores_mean + dem + control + 
                perc_coal + perc_renew + log_cap + log_sales + muni,
              data = subset(agg.df,member_both==0))
summary(reg.3)

reg.4 <- polr(coalition_fac ~ scores_mean + dem + control + perc_coal + perc_renew + log_cap + log_sales + muni 
              + gdp_per2 + per_capita_tax2 + wind_potential_adj + medianaqi + unemp + electricity_price,
              data = subset(agg.df,member_both==0))
summary(reg.4)


#table A6 regressions 
agg.df$coalition_score <- agg.df$num_support - agg.df$num_oppose

reg.1 <- lm(coalition_score ~ scores_mean,
            data = agg.df)
summary(reg.1)

reg.2 <- lm(coalition_score ~ scores_mean + dem + control,
            data = agg.df)
summary(reg.2)

reg.3 <- lm(coalition_score ~ scores_mean + dem + control + 
              perc_coal + perc_renew + log_cap + log_sales + muni,
            data = agg.df)
summary(reg.3)

reg.4 <- lm(coalition_score ~ scores_mean + dem + control + 
              perc_coal + perc_renew + log_cap + log_sales + muni
            + gdp_per2 + per_capita_tax2 + wind_potential_adj + medianaqi + unemp + electricity_price,
            data = agg.df)
summary(reg.4)





